In [1]:
import warnings
warnings.filterwarnings("ignore")

Configuration¶

In [2]:
import scanpy as sc
import scvi
import anndata as ad
import torch
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os

# Make sure plots appear inline in notebooks
%matplotlib inline
# Set Scanpy plotting settings
sc.settings.set_figure_params(dpi=100, facecolor='white', frameon=False)
# Set scvi-tools settings
scvi.settings.seed = 0
print("scvi-tools version:", scvi.__version__)
print("PyTorch version:", torch.__version__)

# --- Configuration ---
N_HVG = 3000                   # Number of highly variable genes to select
N_NEIGHBORS = 15               # Number of neighbors for KNN graph
N_PCS = 30                     # Number of PCs for KNN graph (used before scVI)
LEIDEN_RESOLUTION = 0.6        # Resolution for Leiden clustering
SCVI_MAX_EPOCHS = 400          # Maximum training epochs for scVI (can be lowered for speed)

Seed set to 0

scvi-tools version: 1.3.0
PyTorch version: 2.6.0+cu124
In [3]:
sampleid  = "Explanted1"
results_dir = "./102_Scvi_visium_results_"+sampleid+"/"
os.makedirs(results_dir, exist_ok=True)


adata = sc.read_h5ad("./101_Preprocess_processed_adata/QC_processed_"+sampleid+".h5ad")

print("Original AnnData object:")
print(adata)

Original AnnData object:
AnnData object with n_obs × n_vars = 850 × 18074
    obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'pass_qc'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
    uns: 'spatial'
    obsm: 'spatial'

Step 1: Preserving raw counts...¶

In [4]:
if 'counts' not in adata.layers:
    adata.layers['counts'] = adata.X.copy()
else:
    print("Layer 'counts' already exists. Assuming it contains raw counts.")

Step 2: Basic Quality Control and Filtering¶

In [5]:
# 2a. Filter out spots not 'in_tissue' (usually done by default by spaceranger/read_visium)
if 'in_tissue' in adata.obs.columns:
    n_spots_before = adata.n_obs
    adata = adata[adata.obs['in_tissue'] == 1, :]
    print(f"  Filtered out {n_spots_before - adata.n_obs} spots not in tissue.")
else:
    print("  'in_tissue' column not found, assuming all spots are in tissue.")

  Filtered out 0 spots not in tissue.
In [6]:
adata.layers['counts'] = adata.X.copy() # Ensure layer 'counts' has filtered raw counts
adata.raw = adata.copy() # Store the filtered data with raw counts in .X into .raw

Step 3: QC¶

In [7]:
sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"], jitter=0.4, multi_panel=True)
No description has been provided for this image

Step 4: Prepare Data for scVI¶

In [8]:
# 4a. Normalize total counts per spot and log-transform (for HVG selection and visualization)
adata_for_hvg = adata.copy() # Work on a copy for HVG selection
sc.pp.normalize_total(adata_for_hvg, target_sum=1e4)
sc.pp.log1p(adata_for_hvg)
In [9]:
# 4b. Identify Highly Variable Genes (HVGs) using the log-normalized data
print(f"  Selecting top {N_HVG} Highly Variable Genes...")
sc.pp.highly_variable_genes(
    adata_for_hvg,
    n_top_genes=N_HVG,
    subset=False, # Don't subset yet, just mark them
    flavor="seurat_v3", # Common flavor, works well
    layer=None # Use adata_for_hvg.X which is log-normalized
)

  Selecting top 3000 Highly Variable Genes...
In [10]:
# Transfer the HVG information back to the original adata object
adata.var['highly_variable'] = adata_for_hvg.var['highly_variable']
adata.var['highly_variable_rank'] = adata_for_hvg.var['highly_variable_rank']
adata.var['means'] = adata_for_hvg.var['means']
adata.var['variances'] = adata_for_hvg.var['variances']
adata.var['variances_norm'] = adata_for_hvg.var['variances_norm']

print(f"  Marked {adata.var['highly_variable'].sum()} genes as highly variable.")

  Marked 3000 genes as highly variable.

Step 5: Setup and Train scVI Model¶

In [11]:
# 5a. Setup AnnData object for scvi-tools
# Tell scVI to use the raw counts stored in the 'counts' layer
# It will automatically use the 'highly_variable' column in adata.var if setup correctly
scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts", # Use the raw counts layer we created
    # Optional: Add batch key if needed, e.g., batch_key="sample_id"
    # Optional: Add categorical covariates if relevant, e.g., categorical_covariate_keys=["cell_type"]
)
In [12]:
# 5b. Initialize the scVI model
# n_latent: Dimensionality of the latent space (adjust if needed, 10-30 often works well)
# n_layers: Number of hidden layers in the neural network (1 or 2 is common)
vae = scvi.model.SCVI(adata, n_layers=2, n_latent=30, gene_likelihood="zinb") # ZINB recommended for UMI counts
In [13]:
# 5c. Train the scVI model
# use_gpu=USE_GPU will automatically use GPU if available and specified
# plan_kwargs can control learning rate, weight decay, etc.
# check_val_every_n_epoch=10 helps with early stopping monitoring
print(f"  Training scVI model for up to {SCVI_MAX_EPOCHS} epochs...")
vae.train(max_epochs=SCVI_MAX_EPOCHS, plan_kwargs={'lr': 1e-3}, check_val_every_n_epoch=10)

GPU available: True (cuda), used: True

TPU available: False, using: 0 TPU cores

HPU available: False, using: 0 HPUs

You are using a CUDA device ('NVIDIA GeForce RTX 4060 Ti') that has Tensor Cores. To properly utilize them, you should set `torch.set_float32_matmul_precision('medium' | 'high')` which will trade-off precision for performance. For more details, read https://pytorch.org/docs/stable/generated/torch.set_float32_matmul_precision.html#torch.set_float32_matmul_precision

LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

  Training scVI model for up to 400 epochs...

`Trainer.fit` stopped: `max_epochs=400` reached.
In [14]:
# Extract ELBO loss from training history
history = vae.history
elbo_train = history["elbo_train"]  # Training loss
elbo_validation = history["reconstruction_loss_train"]  # Validation loss

# Plot convergence curve
plt.figure(figsize=(6, 4))
plt.plot(elbo_train, label="Training Loss", color="blue")
plt.plot(elbo_validation, label="Validation Loss", color="red", linestyle="dashed")
plt.xlabel("Epochs")
plt.ylabel("Negative ELBO")
plt.title("scVI Training Convergence")
plt.legend()
plt.savefig(os.path.join(results_dir, "scvi_training_elbo.png"))
No description has been provided for this image

Step 6: Post-Training Analysis - Latent Space, Clustering, Visualization¶

In [15]:
# 6a. Get the latent representation from the trained model
print("  Extracting scVI latent representation...")
adata.obsm["X_scVI"] = vae.get_latent_representation()

  Extracting scVI latent representation...
In [16]:
# 6b. Perform clustering on the scVI latent space
print("  Performing Leiden clustering on scVI latent space...")
sc.pp.neighbors(adata, n_neighbors=N_NEIGHBORS, use_rep="X_scVI")
sc.tl.leiden(adata, resolution=LEIDEN_RESOLUTION, key_added="leiden_scvi", flavor="igraph", n_iterations=2)

  Performing Leiden clustering on scVI latent space...
In [17]:
# 6c. Compute UMAP embedding based on the scVI latent space
print("  Computing UMAP embedding...")
sc.tl.umap(adata, min_dist=0.3) # min_dist controls spread

  Computing UMAP embedding...
In [18]:
# 6d. Visualize the results
print("  Generating visualizations...")

# UMAP colored by cluster
sc.pl.umap(adata, color="leiden_scvi", title="UMAP colored by scVI Leiden Clusters", show=True)

  Generating visualizations...
No description has been provided for this image
In [19]:
# Spatial plot colored by cluster
sc.pl.spatial(adata, color="leiden_scvi", title="Spatial - scVI Leiden Clusters", show=True, spot_size=140) # Adjust spot_size as needed
No description has been provided for this image
In [20]:
# Visualize QC metric spatially (optional)
sc.pl.spatial(adata, color="total_counts", title="Spatial - Total Counts", show=True, spot_size=140, cmap="jet")
No description has been provided for this image
In [21]:
sc.pl.spatial(adata, color="n_genes_by_counts", title="Spatial - Genes per Spot", show=True, spot_size=140)
No description has been provided for this image

Step 7: Differential Gene Expression (DGE) using scVI¶

In [22]:
de_df = vae.differential_expression(groupby="leiden_scvi")
In [23]:
de_df
Out[23]:
proba_de proba_not_de bayes_factor scale1 scale2 pseudocounts delta lfc_mean lfc_median lfc_std ... raw_mean1 raw_mean2 non_zeros_proportion1 non_zeros_proportion2 raw_normalized_mean1 raw_normalized_mean2 is_de_fdr_0.05 comparison group1 group2
PNPLA2 0.9390 0.0610 2.733941 3.085768e-03 8.745411e-04 0.000384 0.25 2.171902 2.321421 1.144848 ... 3.377049 2.851645 0.909836 0.836538 34.246891 8.862825 False 0 vs Rest 0 Rest
PLIN4 0.9378 0.0622 2.713182 2.772105e-03 5.255678e-04 0.000384 0.25 3.554765 3.463848 2.212646 ... 2.836066 0.967035 0.852459 0.440934 30.232544 5.199165 False 0 vs Rest 0 Rest
ADIPOQ 0.9216 0.0784 2.464287 3.732225e-03 6.853951e-04 0.000384 0.25 3.688207 3.605845 2.448200 ... 3.901640 1.101650 0.885246 0.461538 40.987972 6.746245 False 0 vs Rest 0 Rest
ADH1B 0.9190 0.0810 2.428837 1.494951e-03 3.714916e-04 0.000384 0.25 2.555676 2.617056 1.505686 ... 1.434425 1.133245 0.696721 0.579670 16.270020 4.085680 False 0 vs Rest 0 Rest
CIDEC 0.9162 0.0838 2.391802 6.816586e-04 1.209419e-04 0.000384 0.25 3.576041 3.673707 2.115006 ... 0.590164 0.188187 0.508197 0.148352 7.030323 1.026994 False 0 vs Rest 0 Rest
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
HEATR9 0.0000 1.0000 -0.000000 2.022143e-07 1.258345e-06 0.000328 0.25 -0.723028 -0.265741 1.443112 ... 0.005814 0.008850 0.005814 0.008850 0.039578 0.007806 False 6 vs Rest 6 Rest
STEAP1B 0.0000 1.0000 -0.000000 2.536143e-07 1.193080e-06 0.000328 0.25 -0.533144 -0.139834 1.406024 ... 0.000000 0.002950 0.000000 0.002950 0.000000 0.003385 False 6 vs Rest 6 Rest
STK31 0.0000 1.0000 -0.000000 1.673650e-07 1.044207e-06 0.000328 0.25 -0.528764 -0.224700 1.245601 ... 0.000000 0.005900 0.000000 0.005900 0.000000 0.023605 False 6 vs Rest 6 Rest
CCDC114 0.0000 1.0000 -0.000000 2.915561e-07 5.002540e-07 0.000328 0.25 0.148614 0.280663 1.032658 ... 0.005814 0.002950 0.005814 0.002950 0.068886 0.001761 False 6 vs Rest 6 Rest
SLC51A 0.0000 1.0000 -0.000000 8.948022e-07 9.602747e-07 0.000328 0.25 0.101715 0.161527 1.311177 ... 0.000000 0.005900 0.000000 0.005900 0.000000 0.008950 False 6 vs Rest 6 Rest

126518 rows × 22 columns

In [24]:
de_df = de_df[de_df.proba_de > 0.85]
de_df = de_df[de_df.lfc_mean > 1.0]
In [25]:
de_df = de_df.sort_values(by=["group1", "proba_de", "lfc_mean"], ascending=[True, False, False])

dge_filename = os.path.join(results_dir, "scvi_differential_expression.csv")
de_df.to_csv(dge_filename)
print(f"  Differential expression results saved to '{dge_filename}'")

  Differential expression results saved to './102_Scvi_visium_results_Explanted1/scvi_differential_expression.csv'
In [26]:
print("\n  Top markers per cluster (based on scVI DGE):")
n_top_markers = 5
for cluster_id in sorted(de_df['group1'].unique()):
    print(f"\n  --- Cluster {cluster_id} vs Rest ---")
    top_markers = de_df[de_df['group1'] == cluster_id].head(n_top_markers)
    if top_markers.empty:
        print("    No significant markers found with current thresholds.")
    else:
        print(top_markers[['proba_de', 'lfc_mean', 'non_zeros_proportion1']]) # Show key stats

  Top markers per cluster (based on scVI DGE):

  --- Cluster 0 vs Rest ---
        proba_de  lfc_mean  non_zeros_proportion1
PNPLA2    0.9390  2.171902               0.909836
PLIN4     0.9378  3.554765               0.852459
ADIPOQ    0.9216  3.688207               0.885246
ADH1B     0.9190  2.555676               0.696721
CIDEC     0.9162  3.576041               0.508197

  --- Cluster 1 vs Rest ---
       proba_de  lfc_mean  non_zeros_proportion1
RGS5     0.9518  2.110264               1.000000
TPM2     0.9404  2.430894               1.000000
FLNA     0.9398  2.155134               1.000000
MRVI1    0.9396  2.725575               0.994118
ACTN1    0.9356  1.934578               1.000000

  --- Cluster 2 vs Rest ---
        proba_de  lfc_mean  non_zeros_proportion1
BGN       0.9442  2.300829               0.958904
FN1       0.9412  1.968425               0.986301
VCAN      0.9386  2.378482               0.958904
S100A4    0.9382  1.639124               0.849315
COL3A1    0.9264  1.403500               0.904110

  --- Cluster 3 vs Rest ---
       proba_de  lfc_mean  non_zeros_proportion1
FBLN1    0.8506  2.477116               0.908257

  --- Cluster 5 vs Rest ---
       proba_de  lfc_mean  non_zeros_proportion1
CCL14    0.8594  3.204554               0.715789

Step 8: Visualizing top marker genes spatially (optional)...¶

In [27]:
# Get some top markers across a few clusters
markers_to_plot = []
for cluster_id in sorted(de_df['group1'].unique())[:3]: # Plot for first 3 clusters
    cluster_markers = de_df[de_df['group1'] == cluster_id].head(2).index.tolist()
    if cluster_markers:
        markers_to_plot.extend(cluster_markers)

markers_to_plot = list(dict.fromkeys(markers_to_plot)) # Unique markers
In [28]:
markers_to_plot
Out[28]:
['PNPLA2', 'PLIN4', 'RGS5', 'TPM2', 'BGN', 'FN1']
In [29]:
sc.pl.umap(adata, color = markers_to_plot, cmap='jet')
No description has been provided for this image
In [30]:
if markers_to_plot:
    print(f"  Plotting spatial expression for: {markers_to_plot}")
    # Use expression from adata.raw for visualization if desired (raw counts)
    # Or use normalized data from adata.X
    plt.figure()
    sc.pl.spatial(
        adata,
        color=markers_to_plot,
        spot_size=180,
        cmap = 'jet',
        alpha=1,
        ncols=min(len(markers_to_plot), 4), # Adjust layout
        # layer='counts', # Uncomment to plot raw counts from the layer
        use_raw=True,
        show=True,
    )
    plt.close()
    print(f"  Marker gene spatial plots saved with prefix in '{results_dir}/figures/'")
else:
    print("  No top markers found to plot.")

  Plotting spatial expression for: ['PNPLA2', 'PLIN4', 'RGS5', 'TPM2', 'BGN', 'FN1']

<Figure size 400x400 with 0 Axes>
No description has been provided for this image 
  Marker gene spatial plots saved with prefix in './102_Scvi_visium_results_Explanted1//figures/'

Step 9: Saving final AnnData object and scVI model...¶

In [31]:
adata_filename = os.path.join(results_dir, "processed_visium_adata_scvi.h5ad")
adata.write(adata_filename)
print(f"  Processed AnnData object saved to '{adata_filename}'")

model_filename = os.path.join(results_dir, "scvi_model")
vae.save(model_filename, overwrite=True)
print(f"  Trained scVI model saved to '{model_filename}'")

  Processed AnnData object saved to './102_Scvi_visium_results_Explanted1/processed_visium_adata_scvi.h5ad'

  Trained scVI model saved to './102_Scvi_visium_results_Explanted1/scvi_model'